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Abstract 

We consider systems of ordinary differential equations with known first 
integrals. The notion of a discrete tangent space is introduced as the or- 
thogonal complement of an arbitrary set of discrete gradients. Integrators 
which exactly conserve all the first integrals simultaneously are then de- 
fined. In both cases we start from an arbitrary method of a prescribed 
order (say, a Runge-Kutta scheme) and modify it using two approaches: 
one based on projection and one based one local coordinates. The meth- 
ods are tested on the Kepler problem. 



1 Introduction 

A system of ordinary differential equations which preserves a first integral H (y) 
can be written in the form 



V = f(y) = S(y)WH(y), y E R"\ (1.1) 

where S(y) is an antisymmetric matrix, see |12| . An approximate numerical 
solution, y n s» y(t n ), n > 1, is said to be integral preserving if H(y n ) — 
H(y°), n > 1. Energy preserving methods go all the way back to the seminal 
paper by Courant, Friedrichs, Lewy [4], where an energy preserving difference 
scheme is derived and the property is used to prove convergence of the scheme. 
Recently, it has become increasingly popular to study conservative schemes as 
discrete dynamical systems in their own right, attempting to mimic properties 
of the continuous system by introducing suitable discrete counterparts. An ex- 



ample of importance in this note is the replacement of the gradient in (1.1 ) by 
a discrete gradient operator. 

The discrete gradient method is one of the most prevalent approaches in the 
literature. It was first systematically treated by Gonzalez [5] and McLachlan 
et al. [12] , The idea is to introduce a discrete approximation to the gradient, 
letting VH : R m x R m — > W n be a continuous map satisfying 

H(u) - H{v) = VH(v, u) T (u - v), 
VH(u,u) = VH(u). 



The existence of such discrete gradients is well established in the literature, see 
for instance the monograph by Hairer et al. [5] or the papers [S] and |12j . Their 



construction is not unique, we give here two different examples. The Averaged 
Vector Field (AVF) gradient [T7] is defined as 

W AVF H(v,u)= [ VZT(£« + (l-0«)d£. (1.2) 
Jo 

The coordinate increment method [7J is defined in terms of the coordinates of 
the vectors v and u, the ith component of ^7H(v,u) is then given as 



« „, u ^(«l,---,Ui,«i+l,---,Vm)- H(u 1 ,...,v, i _i,v i ,...,v m ) 
{VciH{v,u))i = . 



An important difference between these two discrete gradients is that (1.2) is 



(1.3) 



symmetric, V avfH(v,u) = Vavf#(u, v), while (1.3) is not. However, note 
that a symmetric version of the coordinate increment discrete gradient can be 
constructed by 

V SC iff(v, u ) = \ (Vci-ff(«, u) + VciH (u, «)) . (1.4) 

Once a discrete gradient has been found, one immediately obtains an integral 
preserving method by simply letting 

,,n+l _ ,.n _ 

V h ~ = S(y n ,y n+1 )\7H(y n ,y n+1 ). 

Here h is the time step, and S(y n ,y n+1 ) is some skew-symmetric approximation 
to the matrix S, one would normally require that S(y) = S(y,y). We remark 
that discrete gradient methods are implicit. 

In this note we consider the case where there are more than one first in- 
tegral and the objective is to preserve any number of such invariants simulta- 
neously. Some earlier attempts to achieve this include the papers [12l [18] in 
which the antisymmetric matrix S(y) is replaced by an antisymmetric tensor 
taking discrete gradients of all integrals to be preserved as input. A formula 
for this antisymmetric tensor is given. Another approach is an integrator for 
a class of separable Hamiltonian systems (P3] and references therein), where 
the integrator which preserves all integrals is designed based on separation of 
variables by the Kustaanheimo-Stiefel transformation [20 . This transformation 
was also adopted in the development of the "exact" integrator for the Kepler 
problem by Kozlov ([8], see also [I]). It would be also noteworthy that Labudde 
and Greenspan [9 proposed an energy-and-angular-momentum-preserving inte- 
grator for the differential equations of motion of classical mechanics, and they 
developed similar integrators of high order of convergence in the sequels [10[ 111) . 
Energy and linear momentum are exactly conserved, but for the system case |llj 
angular momentum is not. Another approach is used by Simo et al. in [19 to 
develop schemes that preserve energy and momentum. 

We shall instead present an approach which does not rely on finding such a 
tensor nor a structure of the equation, we only assume knowledge of the first 
integrals to be preserved as well as the ODE vector field itself. Examples of 
such invariants are the energy and momentum, but our approach is not limited 
to these. The discrete gradients are essential to the algorithm we develop, but 
note that the methods themselves are not discrete gradient methods in the 
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usual sense. After defining the general method, we present two particular cases 
based on projection and local coordinates, respectively. Both use an underlying 
scheme of arbitrary order p, we prove that the resulting schemes retain this 
order, regardless of the choice of discrete gradient. 

In Section[3]we apply the new methods to the Kepler problem, a system with 
four degrees of freedom and three independent first integrals. We illustrate our 
approach by preserving combinations of one or more of these three integrals. 
An interesting question is whether the present methods based on discrete gra- 
dients perform better than the more standard projection methods which make 
use of the exact gradients of the first integrals. We show two examples where 
the different approaches are compared. In the examples we use Runge-Kutta 
methods of different order as the underlying schemes. The Kepler problem is of 
course a well-known and popular test case, and methods which preserve one or 
more integrals for this particular problem can be found in e.g. [21 [3l [13], the 
methods in these references are derived by means of the Kustaanheimo-Sticfcl 
transformation. 

2 Preserving multiple invariants 

Suppose that an ODE system ( |1.1[ ) possesses q > 1 independent first integrals, 
Hi(y), . . . , H q (y). These invariants foliate R m into (m — g)-dimensional sub- 
manifolds (leaves) 

M = M c = {ye R m : H^y) = c u H 2 (y) = c 2 , . . . , H q (y) = cj. 

The tangent space T y M of M at y is the orthogonal complement to 

span{Vifi (y),... ,VH q (y)}. 

For simplicity we write only M for M c for the rest of this paper. 

Definition 2.1. Let V be a fixed discrete gradient operator and let Hi, . . . , H q 

be independent first integrals. The discrete tangent space at (it, v) £ K m x R m 
is 

T M M = {neR m : (VHj(v, u), if) =0, 1 < j < q}. 
A vector r\ = rj^ v u ^ e T( V ^M is called a discrete tangent vector. 
Note that this definition causes T( y ^M = T y M. 
Lemma 2.2. Any integrator satisfying 

!/" +1 -!/"^ (f ^)^(f^) M 
preserves all integrals, in the sense that Hi(y n+1 ) — Hi(y n ), 1 < i < q. 
Proof. For any i we compute 

H t (y n+1 ) - H 2 (y n ) = Vi/ 4 (y", y n+1 ) T (y n+1 - y n ) 

= WH l (y n ,y n+1 ) T V(y^ y ^)=0. 

□ 
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This way of devising integral preserving schemes is similar, but slightly dif- 
ferent to the presentation in e.g. |BJ, [H] and [T7]. In this paper we will outline 
two ways of ensuring that the condition of Lemma 2.2 is satisfied - projection 
and local coordinates. We emphasise, however, that there are other ways of 
satsifying Lemma 2.2 One example taken from [T^] is 

= hS(y n lV n+1 WH l {y n ,y n+1 ) . . . VH q (y n , y n+1 ), 
where S is a g-dimensional skew-symmetric tensor. 

2.1 Projection 

We consider with the first integrals H\{y), . . . , H q {y). We propose the 

projection scheme 



,n+i _ 



Mv n ), y 



n+l n 



(2.1) 



where 4>h is the discrete flow that defines an arbitrary method of order p 

y(t + h)- u n+1 = y(t + h)- 4>h{y{t)) = 0(hP +1 ), 

and V(y n , y n+1 ) is a smooth projection operator onto the discrete tangent space 
Ttyn yn+\}M. An alternative method is 



y 



n+i = y n + hv(y n ,y n+1 )My n ,y n+1 ) (2-2) 



» n+1 = y n + hMy n ,y n+1 )- 



where iph is the increment function that defines a method of the form 

y" 

This method is itself assumed to be of order p, that is 

y(t + h) - y(t) - h1> h {y{t),y{t + h)) = 0(h^ +1 ). (2.3) 

We remark in passing that if both the discrete gradients V Hi and the increment 
function iph are symmetric, then the method (2.2) is symmetric. 

Example 2.3. Using Runge-Kutta as the underlying scheme <ph we can con- 
struct examples of (2.1 1. The unprojected solution is given as 

s 

y n + hJ2b t h, 



n+l 



where k\, . . . , k s are the solutions to the (possibly implicit) equations 

( 

\ 3=1 



Even if the Runge-Kutta scheme is explicit the scheme (2.1 1 will be implicit 



since y 



n+l 



appears in the projection operator 



,,n+l 



r + hViy^y^l^hki 
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The difference between (2.1) and (2.2 1 is subtle, but to illustrate that they 
are in fact distinct we consider the implicit midpoint method as the underlying 
scheme and we get for the two methods 

n+1 =y n + hV(y n ,y n+1 )fl y 



II 
II 

where in the former method u n+1 is computed by solving 



' -V n + hV(y n ,y n+1 )fl y 1 



u n +t =y n + hf 



y n + u n+l 



Theorem 2.4. The schemes (2.1 1 and (2.2) are of order p, that is 
y(t + h)- y(t) - V(y(t),y(t + h))(u n+1 - y(t)) = 

u n+1 = Mv(t)l 

and 



(2.4) 



y(t + h)- y(t) - hV(y(t), y{t + h))^ h (y(t), y(t + h)) = (2.5) 

Proof. We use the shorthand notation V for V(y(t),y(t + h)) in this proof. To 
prove (2.4), we compute 

y(t + h) - y{t) - V(u n+1 - y(t)) 
= y(t + h)- y(t) -{!-{!- T))(u n+1 - y(t)) 
= y(t + h) - y(t) - (u n+1 - y{t)) + (1 - V)(u n+1 - y(t)). 



Since 4>h is of order p, we have 



y(t + h) - y(t) - (■ 



,n+l 



/(*)) = y(t + h)- u n+1 = o(h p+1 ). 



(2.6) 



Therefore if we have 



(l-V)(u n+1 -y(t))^0(hP +1 ), 



the proof is completed. This estimate is obtained in the following way. Because 
the image of I — V(y(t), y(t + h)) is spanned by 

{V#i(y(t),y(i + h)), . . .,VH q (y(t),y(t + /»))}, 

it is enough to show 

VHi(y(t),y(t + h)) ■ K+ 1 - y{t)) = 0(h" +1 ). 



From (2.6) we obtain 

u n+1 - y(t) = y(t + h)- y(t) + 0(^ +1 ), 
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and hence 



VF i (y(t),y(t + / l ))-K +1 -t/(t)) 
= VHMt),y(t + h)) ■ (y(t + h)- y(t) + 
= Hi(y(t + h)) - Hi(y(t)) + 
= 0{h p+1 ). 

The last equality is from the conservation property of the original equation. The 



proof of (2.5) is almost identical and therefore omitted. □ 



We remark that in what we call standard projection methods u n+1 is projected 
orthogonally onto the manifold M by computing 

min ||y n+1 - subject to y n+1 £ M. 

This can be achieved for instance by using Lagrange multipliers (see [BJ, for 
instance). This approach differs from ours since we project along the discrete 
gradients which depend on the end point y n+1 . 

Computing the projector. A simple and straightforward way of obtain- 
ing the projector V(y n ,y n+1 ) is as follows: Define the (q x m)-matrix Y — 
Y(y n , y n+1 ) whose columns are the discrete gradients VHi(y n , y n+1 ), i = 1, . . . , q 
Compute a reduced Q ^-decomposition QR = Y where Q E M. rnxq and R £ 
MS*?. Then define the projection matrix as V(y n ,y n+1 ) = I - QQ T . 

2.2 Local coordinates 

The local coordinates approach presented here is basically of the same type 
as the standard method by Potra and Rheinboldt ([IS], see also [SJ HE])- One 
important difference is that our local coordinates are algorithmically constructed 
by using discrete gradients. We also present an "automatic differentiation" 
algorithm of the coordinate map, which could be used to increase the efficiency 
of the computations. 

Inspired by PQ, we consider local coordinates on a chart containing y° by 
defining a map r/ 1— > y = xi 7 ])- The map is defined implicitly by 

X (ri) = y- y-y° = T(y° 1 y)r ]l (2.7) 

where T(y°,y) is a smooth m x (m — g)-matrix whose columns form a basis 
for the left nullspace (orthogonal column complement) of the matrix Y(y°,y) = 
\VHi(y°,y), . . . ,VH q (y° ,y)]. We suppress the dependency on y° and use the 
shorthand notation T(y) and Y(y) for the rest of this paper. 

Lemma 2.5. Suppose that Vi?i(y ), . . . , \?H q (y°) are linearly independent for 
all y° e M. Suppose also that for all y° G M, VHi(y°, y), . . .,VH q (y°, y) are 
C°° with respect to y. Then the following statements hold. 



1. (2.7) defines a one-to-one map uj y o in a neighborhood N y o C M. ui y a ana 
„ oreC™. 



CJ o 
V 



2. The collection of the pairs {(N y o,oj y o) | y° € M} forms an atlas of M. 
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Proof. 1. From the continuity of the discrete gradients, we deduce that for 
all y° £ M, there exists a neighborhood N y o in M. m of y° in which 
VHi(y°,y), . . . ,VH q (y°,y) are linearly independent. For y £ N y o, T(y) 
admits the QR decomposition T(y) — QR and rj is obtained by 77 = 
(R T R)- 1 Q T (y-y°) = (T T T)~ 1 r T (y — y°). This is a C°° function. Con- 
versely, the Jacobian matrix of the function r](y) at y = y° is 

fV) = (T T T) _1 T T 
ay 



and hence 

rfln\ 

<9y 



rank|^(y°) = rank(T T T) -:l T T = m - q. 



Thus to is defined in a neighborhood N y o of y° by the implicit function 
theorem and is C°° . The proof is completed by letting N y a = N y o n N y o n 
M. 

2. This is immediately obtained from the first statement. 

□ 

Consider now the curve r/(t) and let y(t) = x(ry(t)). We differentiate the 
curve to obtain from (2.7) 

y(t) = T^ t) (y(t))r 1 (t)+T(y(t))r 1 (t). 

From this we compute 

V = -T T ( X o vK ov (f(x o r)))V + T T ( X o V )f(x ° V) (2-8) 

where the original ODE is y = f{y). The method we propose takes one step as 
follows 

1. Let rjo = 0. 



2. Take a step with any pth order method applied the ODE (2.8| using 
y° = y n in p?7| . The result is r? x . 

3. Compute = x( r ?i)- 



We immediately obtain the next theorem from Lemma 2.5 because the solution 
curve lies in M and a pth order method is applied in a chart of M. 



Theorem 2.6. Under the assumptions of Lemma 2.5 the above scheme is of 
order p. 

The main difficulty in this approach is the computation of the derivative map 
T y {C) for arbitrary values of y £ W n and ( £ ]R m . This is needed explicitly in the 
integration algorithm, but may also be a useful tool in computing the coordinate 
map ( |2.7[ ). We may define T(y) as the last m— q columns of the m x m-matrix 
Q{y) defined through a QR-decomposition where Y(y) = Q(y)R(y) and where 
we have used the shorthand notation 

Y(y) = [\7H 1 (y n ,y),...,\7H q (y n ,y)} 
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We realise the QR-decomposition by means of the Householder method, 
applying a sequence of q elementary orthogonal transformations to the matrix 
Y(y) as described in most elementary text books in numerical linear algebra, 
see e.g [21]. Each transformation is of the form 

Q k =I- 2v k vJ, v k £ R m , vjv k = 1, 

and its aim is to eliminate all elements under the diagonal in the fcth column of 
the matrix to which it is applied. 

In order to explain how we compute the derivative Q' y (C) =: DQ(y,£), we 
first review the Householder method. 

y(i) .-y 
for k = 1 : q, 



r«-||n fe >f>||e fc 



Vk - — - 



for r = k : q, 

Y r (k+1) = (I-2v k vJ)Y r {k) 
end 
end 

where the following conventions have been used: 

• || • || is the Euclidean norm. 

• F( fe+1 ) = Q fc y( fc \ y r (fe) is column r of Y^ . 

• The projector puts zeros in the first k — 1 components and leaves the 
rest of the components unchanged when applied to a vector in R m . 

• e k is the fcth canonical unit vector in K m . 

The vectors v k computed in the algorithm contain all information needed 
to reconstruct the factor Q, whereas R :— Y^- q+1 \ For simplicity, and to avoid 
the loss of regularity in Q(y) viewed as a matrix valued function of y, we have 
here ignored the sign convention which is usually applied in the definition of w k 
[2"T] . The idea is now to differentiate the variables in the algorithm with respect 
to y, writing for any object, say X(y), its derivative as DX — DX(y,() — 
^\ s _ X(y + e() for any y, ( g E m . The dependence on y, ( will usually be 
suppressed in the notation when no confusion is at risk. Notice that D commutes 
with Tl k for any fc. The following recursion formulae are easily derived 

Dw — n DY« (UkYtYUkDY^ 
L>w k — n k L>r k t— e fc , 

lin^H 

_. wjDw k \ 

Dv k = Dw k - ., w k \ w k \ 

V If* II / 

D y r (fe+i) = DY (k) _ 2 (yT Y (k) Dvk + DvjY^v k + vjDY^v^ . 

The initial DY^ = DY^(y, C) should be computed by differentiating the q 
discrete gradients Viii (p, y), . . . , \7H q (p, y) with respect to y. When the discrete 
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gradients are given by the AVF formula, we may derive the following expressions 
for the rth column of DY^ 

DY r W(yX) = DVH r (y,0= (7 V 2 #r(£y + (1 - £>)d^) • C 

where V 2 H r is the Hessian of the integral H r . The expression for the coordinate 
increment cases are given in the appendix. 

In the present application we only make use of the Q-part of the QR- 
decomposition, thus we need only store v% , . . . , v q and Dv\ , . . . , Dv q for sub- 
sequent use. We may summarize the extended algorithm for computing these 
quantities as follows, using a Matlab inspired indexing notation where the sub- 
matrix Y a -j, c:( i of Y means 



Y. 



a:b,c:d 



(Y a , 
\ Y b,, 



Y, 



Y 



b.di 



Given Y and DY as m x q-matrices 
for k = 1 : q 

W = Yk;rn,k ~ \\Yk:m,k^l 

r, ^ v Y k k DYk: m ,k 

Dw = DY k . mk — ex 

\\ik:m,k\\ 



Vk 



Dv k = (Dw 



w 1 Dw 



DYk:m,k+l:q — DYk; m ,k+l:q — 2v k V k DYk :m ,k+l:q ~ 2Dv k V k Y k ;m,k+l:q 
~2v k Dv k r Yk- m ,k+l:q 

end 

The first k — 1 entries of the vectors v k , Dv k are zeros, and the remaining 
m — k + 1 entries are contained in v k , Dv k on exit. The complexity of this 
algorithm is 0(mq 2 + q 3 ). 

One should also note that we always multiply Q (DQ resp) by vectors fj 6 
M. m whose first q columns are zero, this may be taken advantage of in the 
implementation. The procedure for computing Qfj by means of v\, . . . , V k is 
described in [3TJ Algorithm 10.3], the cost is 0(mq). We may extend this 
algorithm so that it computes also DQuj given Dv\ , . . . , Dv q . Defining the 
m x TO-matrix P k = Q k ■ ■ ■ Q q , k = 1, . . . , q, we get the downwards recursion 



ffc-i — Qk-iPk, 

and differentiation yields 



Pq — Qqi Pi — Q — Qi 



;</• 



DP, 



k-l 



DQk^Pk + Q k -iDP k 



= -2{Dv k -i vl_ x + v k -i DvJ^Pk + 2v k - 1 vJ_ 1 )DP k 
The following algorithm results for computing cf> = DQ u> 
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ip — u), = 
for k = q : — 1 : 1 

(/) = (/)- 2(vjip Dv k + DvJipVk + vjcfivk) 
ip = ip — 2vjipv k 
end 

The complexity of this algorithm is 0(mq). 



3 Numerical integration of the Kepler problem 

The Kepler two-body problem describes the motion of two bodies which attract 
each other. By placing the first body in the origin, the position (j/i, 2/2) and the 
velocity (2/3,2/4) of the other body are given by the following four-dimensional 
ODE 

Vi = 2/3, 
2/2 = 2/4, 

2/3 = - 7 2 yi 2 \3/2 > 
2/2 

2/4 



(2/?+2/l) 3 / 2 ' 

This system preserves the Hamiltonian 

Hi(y) = - (y! + yl)--== , 
2 Vvi + 2/2 

the angular momentum 

H 2 (y) = 2/12/4 - 2/22/3, 
and the Runge-Lenz-Pauli vector 

H 3 (y) = 2/22/3 - 2/12/32/4 ''" 



2 ; 



TT ( \ 2 2/1 

#4(2/) = 2/12/4 ~ 2/22/32/4 7=,= — 7 

V2/I + 2/2 



Since q = m = 4 any subset of three out of the four invariants is dependent. 
We want to compare schemes that preserve none, one, two, and all of the invari- 



ants above. We use the projection method (2.1 ) with the standard fourth order 



explicit Runge-Kutta method as the underlying scheme. The discrete gradients 



are calculated using ( 1.4 1. The schemes that preserve one of Hi , H3 are denoted 
as RK4Projl and RK4Proj3, respectively. The scheme that preserves both Hi 
and H3 is called RK4Projl3. The original Runge-Kutta scheme preserves nei- 
ther and is denoted RK4. The RK123Proj scheme is omitted from the plot since 
it produces exactly the ellipsis of the Kepler problem. 



The resulting plots from these methods in Figure 3.1 are arranged according 
to the table 
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RK4 


RK4Projl 


RK4Proj3 


RK4Projl3 



Table 3.1: The location of the schemes in the plots of Figure 



34 



The initial values are taken from section 1.2.3 of [5J, 



i/? = l-e, y° 2 =0, y° 3 =0, " 



1-e' 

where the eccentricity is e = 0.6 and the exact solution has period 2tt. The time 
step is h = 0.2 and we integrate for 50000 steps. 



Figure 3.1 shows the numerical solutions. RK4 spirals inwards until it even- 
tually blows up. RK4Projl has a counterclockwise precession effect. The Runge- 
Lenz-Pauli vector has to do with the orientation of the ellipse and RK4Proj3 
does therefore not exhibit this effect, it will however converge to a smaller circle 
around the origin. The solution of RK4Projl3 shows an improvement compared 
to RK4Projl and RK4Proj3. 

This example illustrates that there are cases where the preservation of one or 
more invariants are important to get a numerical solution with good long term 
properties. Not surprisingly, one observes a gradual improvement in the quality 
of the solution as the number of preserved first integrals increases. The extra 
computational effort needed to preserve multiple integrals compared to one is 
almost negligible, in this example the computation took less than 10% longer. 



In Figure 3.2 we plot the global error of four schemes (RK2Projl23, RK4Projl23, 
RK5Projl23, and RK7Projl23) that preserve the four invariants. The under- 
lying schemes are four RK-schemes of order 2, 4, 5, and 7. We see that the 
schemes attain the order of the underlying scheme, which is what we proved in 
Theorem EU 



Figure 3.3 shows a comparison between a standard orthogonal projection 



method as defined at the end of Section 2.1 and the projection method (2.1 1 
Both preserve H\ and Hi simultaneously and use the implicit midpoint method 
as the underlying scheme. Figure |3.3a shows that our projection method is 



more accurate for e = 0.6 while Figure 3.3b shows that the standard projection 



method is more accurate for e = 0.7. In the standard projection method one 
has to compute the distance to the underlying manifold, usually denoted by 
the Lagrange multiplier A, however for our projection scheme this is already 
(implicitly) known. The resulting nonlinear system will have dimension m + q 
(see section IV. 4 in [5]) compared to m for our proposed method. Our imple- 
mentation uses the same nonlinear solver (Matlab's fsolve) for both methods. 
In Figure [3. 3c| we have adjusted h such that both schemes have the same com- 
putation time, in which case one sees that the new method is slightly better 
than the standard projection method even for e = 0.7. 

Conclusion and further work. We have presented a new methodology for 
preserving multiple first integrals in systems of ordinary differential equations, 
using discrete gradients as the underlying tool. By using the notion of a discrete 
tangent space, two methodologies for designing numerical schemes are easily de- 
rived, projection and local coordinates. The resulting algorithms are relatively 
inexpensive compared to well-known algorithms preserving precisely one first 
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0.5 



-0.5 ■ 




-1 




Figure 3.1: The numerical solution (thin line) of the Kepler problem (3.1 ) using 
the schemes of Table [3~T] with h = 0.2. The first 500 steps are shown. The exact 
solution (thick line) is an ellipse with eccentricity e = 0.6. 
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Figure 3.2: The global error of the four schemes RK2Projl23, RK4Projl23, 
RK5Projl23, and RK7Projl23. The dotted lines are reference lines of exact 
order. 



13 




2 
1.5 
1 

0.5 


-0.5 

-1 

-1.5 

-2 
-2 



2 
1.5 
1 

0.5 


-0.5 
-1 
-1.5 



-2 
-2 






yi 



(b) e = 0.7 and h 









(c) e = 0.7 and h = 0.075 (left) /i = 0.05 (right). 



Figure 3.3: The numerical solution (thin line) of the Kepler problem (3.1 ) using 
the midpoint method with standard orthogonal projection (left) and (2.1) with 
the midpoint method as the underlying scheme. The exact solution (thick line) 
is an ellipse with eccentricity e. Both schemes preserves Hi and Hi exactly. 
In the last row (c) the step size h is adjusted such that both schemes have the 
same computation time. ^ 



integral. There are of course several other well known methods for preserving 
multiple invariants, but we believe that the new methods which are based on 
discrete gradients rather than exact ones may be attractive for certain prob- 
lems. Symmetric schemes are easily constructed, no Lagrange multipliers are 
required, and the schemes can be seen as natural generalizations of the popular 
discrete gradient methods. Although the present paper considers only systems 
of ordinary differential equations, the approach taken may be easily adapted to 
partial differential equations to be considered in future work. 
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A Derivative of discrete gradients of the Coor- 
dinate Increment type 

We here give the expression for the derivative map of the discrete gradients 
defined in terms of the coordinate increment method of Itoh and Abe [7] . We 
write (1.3 1 in compact notation as 



(Vci#M)i 



H(u\iv) - HMi-iu) 



1 < i < m 



The jacobian of this map with respect to its second argument is then simply 
the lower-triangular matrix with elements 



^(H(u\iv) - H(u\i- lV )) 



(DVciH(v,u))i 



H(u\ iV ) - ff(u|i_i«) 



{Ui - Vi) 2 







J < « 

3 = * 

i > i 



Similarly, we can compute the Jacobian of the symmetrised version as follows 



(DV S ciH(v,u))i 



2 Ui-Vi 2 (ui-Vi) 2 ■ 

^{Hjv^u) - H{v\ lU )) 

2 ' 

Ui-Vi 
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